Geocomputation with R

R
Creative & Experimental geography
Author

Sungil Park

Published

March 22, 2023

Modified

April 22, 2023

#|warning: False

Packages

library(sf)
Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(raster)
Loading required package: sp
library(spData)
library(spDataLarge)
library(tmap)
library(mapview)
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ dplyr::select()  masks raster::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(grid)
library(gifski)

벡터 데이터

R class의 벡터 데이터와 다른 의미의, 공간 위치데이터의 점,선,면을 나타내는 데이터

#vignette(package = "sf")

World dataset from spData Package

world %>% head()
Simple feature collection with 6 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
Geodetic CRS:  WGS 84
# A tibble: 6 × 11
  iso_a2 name_long conti…¹ regio…² subre…³ type  area_…⁴     pop lifeExp gdpPe…⁵
  <chr>  <chr>     <chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <dbl>   <dbl>
1 FJ     Fiji      Oceania Oceania Melane… Sove…  1.93e4  8.86e5    70.0   8222.
2 TZ     Tanzania  Africa  Africa  Easter… Sove…  9.33e5  5.22e7    64.2   2402.
3 EH     Western … Africa  Africa  Northe… Inde…  9.63e4 NA         NA       NA 
4 CA     Canada    North … Americ… Northe… Sove…  1.00e7  3.55e7    82.0  43079.
5 US     United S… North … Americ… Northe… Coun…  9.51e6  3.19e8    78.8  51922.
6 KZ     Kazakhst… Asia    Asia    Centra… Sove…  2.73e6  1.73e7    71.6  23587.
# … with 1 more variable: geom <MULTIPOLYGON [°]>, and abbreviated variable
#   names ¹​continent, ²​region_un, ³​subregion, ⁴​area_km2, ⁵​gdpPercap
names(world)
 [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
 [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"     
plot(world)
Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
all

sp 데이터는 st_as_sf()로 sf 형식으로 변환

library(sp)

world_sp = as(world, Class = "Spatial")

world_sf = st_as_sf(world_sp)

Plot map

plot(world["continent"])

st_union() 공간정보 합치기

world_asia = world[world$continent == "Asia", ]
asia = st_union(world_asia)

plot(world["pop"], reset = FALSE)
plot(asia, add = TRUE, col = "red")

st_centroid()공간정보의 중심점 계산

cex : symbol size = sqrt(인구수)/1000

plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000  
world_cents = st_centroid(world, of_largest = TRUE)
Warning: st_centroid assumes attributes are constant over geometries
plot(st_geometry(world_cents), add = TRUE, cex = cex)

Highlight

india = world[world$name_long == "India", ]
plot(st_geometry(india), expandBB = c(0.1, 0.1, 0.1, 0.1), col = "gray", lwd = 3)
plot(world_asia[0], add = TRUE)

Geometry types

st_point(c(5, 2, 3, 1)) %>% plot()

multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
st_multipoint(multipoint_matrix) %>% plot()

linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix) %>% plot()

polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
st_polygon(polygon_list) %>% plot()

polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list) %>% plot()

multilinestring_list = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                            rbind(c(1, 2), c(2, 4)))
st_multilinestring((multilinestring_list)) %>% plot()

multipolygon_list = list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
                         list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
st_multipolygon(multipolygon_list) %>% plot()

gemetrycollection_list = list(st_multipoint(multipoint_matrix),
                              st_linestring(linestring_matrix))
st_geometrycollection(gemetrycollection_list) %>% plot()

Simple feature columns (sfc)

st_sfc() 지리 특성을 하나의 컬럼 객체로 합침

point

point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc %>% plot()

polygon

polygon_list1 = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
polygon1 = st_polygon(polygon_list1)
polygon_list2 = list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2)))
polygon2 = st_polygon(polygon_list2)
polygon_sfc = st_sfc(polygon1, polygon2)
polygon_sfc %>% plot()

multilinestring

multilinestring_list1 = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                             rbind(c(1, 2), c(2, 4)))
multilinestring1 = st_multilinestring((multilinestring_list1))
multilinestring_list2 = list(rbind(c(2, 9), c(7, 9), c(5, 6), c(4, 7), c(2, 7)), 
                             rbind(c(1, 7), c(3, 8)))
multilinestring2 = st_multilinestring((multilinestring_list2))
multilinestring_sfc = st_sfc(multilinestring1, multilinestring2)
multilinestring_sfc %>% plot()

geometry

point_multilinestring_sfc = st_sfc(point1, multilinestring1)
point_multilinestring_sfc %>% plot()

sfc 객체는 CRS에 대한 정보를 추가로 저장할 수 있음

CRS(coordinate reference systems, 좌표계시스템)

EPSG : European Petroleum Survey Group, 지도 투영과 datums 에 대한 좌표계 정보 데이터베이스를 제공

st_crs(points_sfc)
Coordinate Reference System: NA
points_sfc_wgs = st_sfc(point1, point2, crs = 4326)
st_crs(points_sfc_wgs)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

위치데이터 + 속성데이터

st_sf() 를 이용하여 sfc와 class sf의 객체들을 하나로 통합할 수 있음

lnd_point = st_point(c(0.1, 51.5))                 # sfg object
lnd_geom = st_sfc(lnd_point, crs = 4326)           # sfc object
lnd_attrib = data.frame(                           # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
)
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object
lnd_sf
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
Geodetic CRS:  WGS 84
    name temperature       date         geometry
1 London          25 2017-06-21 POINT (0.1 51.5)

Raster 데이터

#install.packages("rgdal") 
library(rgdal)
Please note that rgdal will be retired during 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-6, (SVN revision 1201)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
Path to GDAL shared files: C:/Users/sungi/AppData/Local/R/win-library/4.2/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
Path to PROJ shared files: C:/Users/sungi/AppData/Local/R/win-library/4.2/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-0
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") 
new_raster = raster(raster_filepath)

new_raster
class      : RasterLayer 
dimensions : 457, 465, 212505  (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333  (x, y)
extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : srtm.tif 
names      : srtm 
values     : 1024, 2892  (min, max)
plot(new_raster)

RasterLayer class

한개의 층

Default CRS = WGS84 (resolution scale = degrees)

8*8 Raster data

my_raster = raster(nrows = 8, ncols = 8, res = 0.5, 
                   xmn = -2.0, xmx = 2.0, ymn = -2.0, ymx = 2.0, vals = 1:64)
nlayers(my_raster)
[1] 1
## plotting 
plot(my_raster)

RasterBrick class

muliple layers

단일 다중 스펙트럼 위성 파일 (a single multispectral satellite file)

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
r_brick = brick(multi_raster_file)

nlayers(r_brick)
[1] 4
plot(r_brick)

RasterStack class

multiple layer

메모리의 단일 다층 객체 (a single multilayer object in memory)

raster_on_disk = raster(r_brick, layer = 1)
raster_in_memory = raster(xmn = 301905, xmx = 335745,
                          ymn = 4111245, ymx = 4154085, 
                          res = 30)
values(raster_in_memory) = sample(seq_len(ncell(raster_in_memory)))
crs(raster_in_memory) = crs(raster_on_disk) #같은 좌표 입력

r_stack = stack(raster_in_memory, raster_on_disk)
r_stack
class      : RasterStack 
dimensions : 1428, 1128, 1610784, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
names      :   layer, landsat_1 
min values :       1,      7550 
max values : 1610784,     19071 
plot(r_stack)

CRS

지리좌표계

  • 위,경도

  • 각도로거리 측정

투영좌표계

  • “평평한 표면”위의 데카르트 좌표 기반

  • 원점, x,y축

  • m와 같은 선형 측정 단위

CRS in R

  • epsg코드
    • 일반적으로 짧음
    • 잘 정의된 좌표 시스템 하나만을 참조
  • proj4string 정의
    • 투영 우형, 데이텀 및 타원체와 같은 매개변수를 지정할때 더 많은 유연성
    • 다양한 투영 우형 지정, 기존 유형 수정 가능

st_crs()로 좌표계 조회

st_set_crs()로 좌표계 변경

vector_filepath = system.file("vector/zion.gpkg", package = "spDataLarge")
new_vector = st_read(vector_filepath)
Reading layer `zion' from data source 
  `C:\Users\sungi\AppData\Local\R\win-library\4.2\spDataLarge\vector\zion.gpkg' 
  using driver `GPKG'
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
Projected CRS: UTM Zone 12, Northern Hemisphere
## st_read() : read vector dataset in R sf package

st_crs(new_vector)
Coordinate Reference System:
  User input: UTM Zone 12, Northern Hemisphere 
  wkt:
BOUNDCRS[
    SOURCECRS[
        PROJCRS["UTM Zone 12, Northern Hemisphere",
            BASEGEOGCRS["GRS 1980(IUGG, 1980)",
                DATUM["unknown",
                    ELLIPSOID["GRS80",6378137,298.257222101,
                        LENGTHUNIT["metre",1,
                            ID["EPSG",9001]]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]]],
            CONVERSION["UTM zone 12N",
                METHOD["Transverse Mercator",
                    ID["EPSG",9807]],
                PARAMETER["Latitude of natural origin",0,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",-111,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",0.9996,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["False easting",500000,
                    LENGTHUNIT["Meter",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",0,
                    LENGTHUNIT["Meter",1],
                    ID["EPSG",8807]],
                ID["EPSG",16012]],
            CS[Cartesian,2],
                AXIS["(E)",east,
                    ORDER[1],
                    LENGTHUNIT["Meter",1]],
                AXIS["(N)",north,
                    ORDER[2],
                    LENGTHUNIT["Meter",1]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from GRS 1980(IUGG, 1980) to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]

Raster 데이터에서 좌표계

projection() 함수로 확인하거나 설정

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") 
new_raster = raster(raster_filepath) 
projection(new_raster)
[1] "+proj=longlat +datum=WGS84 +no_defs"

변경

new_raster3  <-  new_raster
projection(new_raster3) <-  "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 
                            +units=m +no_defs"

단위

sf 객체 내에 단위가 들어가있음

st_area() [m^2] 단위 같이 반환

set_units(st_object,units) 로 반환 단위 설정

names(world)
 [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
 [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"     
luxembourg = world[world$name_long == "Luxembourg", ]
south_korea = world[world$name_long == "Republic of Korea", ] 


st_area(luxembourg)
2408817306 [m^2]
st_area(south_korea)
99020196082 [m^2]

Raster는 단위에 대한 속성 정보가 없음

벡터 속성 조작

library(tidyr)
library(dplyr)

st_drop_geometry() sf객체에서 속성 데이터만 가져오기

dim(world)
[1] 177  11
world_df = st_drop_geometry(world)
class(world_df)
[1] "tbl_df"     "tbl"        "data.frame"
dim(world_df)
[1] 177  10

Base R, dplyr 구문으로 조작

world %>% 
  filter(continent == "Asia") %>% 
  select(name_long) %>% 
  plot()

world_agg1  <-  aggregate(pop ~ continent, FUN = sum, data = world, na.rm = TRUE)
world_agg1
      continent        pop
1        Africa 1154946633
2          Asia 4311408059
3        Europe  669036256
4 North America  565028684
5       Oceania   37757833
6 South America  412060811
str(world_agg1)
'data.frame':   6 obs. of  2 variables:
 $ continent: chr  "Africa" "Asia" "Europe" "North America" ...
 $ pop      : num  1.15e+09 4.31e+09 6.69e+08 5.65e+08 3.78e+07 ...
class(world_agg1)
[1] "data.frame"

aggregate()함수를 사용해서 속성 데이터를 그룹별로 집계

  • world$pop 은 sf객체가 아닌 숫자형 객체, geometry 정보는 없음
world_agg2  <-  aggregate(world["pop"], by = list(world$continent),
                       FUN = sum, na.rm = TRUE)
class(world_agg2)
[1] "sf"         "data.frame"
class(world['pop'])
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
class(world$pop)
[1] "numeric"

group_by() 사용

world_agg3 <- world %>% 
  group_by(continent) %>% 
  summarize(pop = sum(pop,na.rm = TRUE),n_countries = n())

world_agg3
Simple feature collection with 8 features and 3 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
Geodetic CRS:  WGS 84
# A tibble: 8 × 4
  continent                      pop n_countries                            geom
  <chr>                        <dbl>       <int>                  <GEOMETRY [°]>
1 Africa                  1154946633          51 MULTIPOLYGON (((36.86623 22, 3…
2 Antarctica                       0           1 MULTIPOLYGON (((-180 -89.9, 18…
3 Asia                    4311408059          47 MULTIPOLYGON (((36.14976 35.82…
4 Europe                   669036256          39 MULTIPOLYGON (((26.29 35.29999…
5 North America            565028684          18 MULTIPOLYGON (((-82.26815 23.1…
6 Oceania                   37757833           7 MULTIPOLYGON (((166.7932 -15.6…
7 Seven seas (open ocean)          0           1 POLYGON ((68.935 -48.625, 68.8…
8 South America            412060811          13 MULTIPOLYGON (((-66.95992 -54.…

인구 많은 3개 대륙

world %>% 
  group_by(continent) %>% 
  summarize(pop = sum(pop, na.rm = TRUE), n_countries = n()) %>% 
  top_n(n = 3, wt = pop) %>%
  arrange(desc(pop)) %>%
  plot()

join 가능

world_coffee <-  left_join(world, coffee_data)
Joining with `by = join_by(name_long)`
class(world_coffee)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
names(world_coffee)
 [1] "iso_a2"                 "name_long"              "continent"             
 [4] "region_un"              "subregion"              "type"                  
 [7] "area_km2"               "pop"                    "lifeExp"               
[10] "gdpPercap"              "geom"                   "coffee_production_2016"
[13] "coffee_production_2017"
plot(world_coffee["coffee_production_2017"])

setdiff()로 일치하지 않는 열 식별

setdiff(coffee_data$name_long, world$name_long)
[1] "Congo, Dem. Rep. of" "Others"

stringr 패키지의 str_subset()

library(stringr)
str_subset(world$name_long, "Dem*.+Congo")
[1] "Democratic Republic of the Congo"
coffee_data$name_long[grepl("Congo,", coffee_data$name_long)]  <-  
  str_subset(world$name_long, "Dem*.+Congo")

world_coffee_match <- inner_join(world, coffee_data)
Joining with `by = join_by(name_long)`
nrow(world_coffee_match)
[1] 46

새로운 열 만들기

base R, mutate(), transmute() 구문 사용

world %>% transmute(pop_dens = pop / area_km2) %>% head()
Simple feature collection with 6 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
Geodetic CRS:  WGS 84
# A tibble: 6 × 2
  pop_dens                                                                  geom
     <dbl>                                                    <MULTIPOLYGON [°]>
1    45.9  (((-180 -16.55522, -179.9174 -16.50178, -179.7933 -16.02088, -180 -1…
2    56.0  (((33.90371 -0.95, 31.86617 -1.02736, 30.76986 -1.01455, 30.4191 -1.…
3    NA    (((-8.66559 27.65643, -8.817828 27.65643, -8.794884 27.1207, -9.4130…
4     3.54 (((-132.71 54.04001, -133.18 54.16998, -133.2397 53.85108, -133.0546…
5    33.5  (((-171.7317 63.78252, -171.7911 63.40585, -171.5531 63.31779, -170.…
6     6.33 (((87.35997 49.21498, 86.82936 49.82667, 85.54127 49.69286, 85.11556…

unite() 사용 열 합치기

seperate() 사용 열 분리

world %>%
  unite("con_reg", continent:region_un, sep = ":", remove = TRUE) %>% head()
Simple feature collection with 6 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
Geodetic CRS:  WGS 84
# A tibble: 6 × 10
  iso_a2 name_long      con_reg    subre…¹ type  area_…²     pop lifeExp gdpPe…³
  <chr>  <chr>          <chr>      <chr>   <chr>   <dbl>   <dbl>   <dbl>   <dbl>
1 FJ     Fiji           Oceania:O… Melane… Sove…  1.93e4  8.86e5    70.0   8222.
2 TZ     Tanzania       Africa:Af… Easter… Sove…  9.33e5  5.22e7    64.2   2402.
3 EH     Western Sahara Africa:Af… Northe… Inde…  9.63e4 NA         NA       NA 
4 CA     Canada         North Ame… Northe… Sove…  1.00e7  3.55e7    82.0  43079.
5 US     United States  North Ame… Northe… Coun…  9.51e6  3.19e8    78.8  51922.
6 KZ     Kazakhstan     Asia:Asia  Centra… Sove…  2.73e6  1.73e7    71.6  23587.
# … with 1 more variable: geom <MULTIPOLYGON [°]>, and abbreviated variable
#   names ¹​subregion, ²​area_km2, ³​gdpPercap
world_unite <- world %>%
  unite("con_reg", continent:region_un, sep = ":", remove = FALSE)
world_unite %>% 
  separate("con_reg", c("continent", "region_un"), sep = ":") %>% head()
Simple feature collection with 6 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
Geodetic CRS:  WGS 84
# A tibble: 6 × 11
  iso_a2 name_long conti…¹ regio…² subre…³ type  area_…⁴     pop lifeExp gdpPe…⁵
  <chr>  <chr>     <chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <dbl>   <dbl>
1 FJ     Fiji      Oceania Oceania Melane… Sove…  1.93e4  8.86e5    70.0   8222.
2 TZ     Tanzania  Africa  Africa  Easter… Sove…  9.33e5  5.22e7    64.2   2402.
3 EH     Western … Africa  Africa  Northe… Inde…  9.63e4 NA         NA       NA 
4 CA     Canada    North … Americ… Northe… Sove…  1.00e7  3.55e7    82.0  43079.
5 US     United S… North … Americ… Northe… Coun…  9.51e6  3.19e8    78.8  51922.
6 KZ     Kazakhst… Asia    Asia    Centra… Sove…  2.73e6  1.73e7    71.6  23587.
# … with 1 more variable: geom <MULTIPOLYGON [°]>, and abbreviated variable
#   names ¹​continent, ²​region_un, ³​subregion, ⁴​area_km2, ⁵​gdpPercap

rename(), setNames() 가능

Raster 속성 조작

R의 레스터 객체(raster objects)는 데이터 속성으로 숫자형(numeric), 정수형(integer), 논리형(logical), 요인형(factor) 데이터 유형을 지원하며, 문자형(character)은 지원하지 않음

문자형으로 이루어진 범주형 변수 값(categorical variables’ values)을 가지고 레스터 객체의 속성을 만들고 싶으면

  • 먼저 문자형을 요인형으로 변환 (또는 논리형으로 변환)하고

  • 요인형 값을 속성 값으로 해서 레스터 객체를 만듬

elev <- raster(nrows = 6, # integer > 0. Number of rows 
               ncols = 6, # integer > 0. Number of columns 
               #res = 0.5, # numeric vector of length 1 or 2 to set the resolution 
               xmn = -1.5, # minimum x coordinate (left border) 
               xmx = 1.5, # maximum x coordinate (right border) 
               ymn = -1.5, # minimum y coordinate (bottom border) 
               ymx = 1.5, # maximum y coordinate (top border) 
               vals = 1:36) # values for the new RasterLayer 

elev
class      : RasterLayer 
dimensions : 6, 6, 36  (nrow, ncol, ncell)
resolution : 0.5, 0.5  (x, y)
extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 1, 36  (min, max)
plot(elev, main = 'raster datasets with numeric valeus')

factor 형식으로 변환

grain_order <- c("clay", "silt", "sand")
grain_char <- sample(grain_order, 36, replace = TRUE)
grain_fact <- factor(grain_char, levels = grain_order)
grain <- raster(nrows = 6, ncols = 6, res = 0.5, 
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
               vals = grain_fact)
plot(grain)

factor 추가

levels(grain)[[1]] <- cbind(levels(grain)[[1]], wetness = c("wet", "moist", "dry"))
levels(grain)
[[1]]
  ID VALUE wetness
1  1  clay     wet
2  2  silt   moist
3  3  sand     dry
grain[c(1, 11, 35)]
[1] 3 3 3
factorValues(grain, grain[c(1, 11, 35)])
  VALUE wetness
1  sand     dry
2  sand     dry
3  sand     dry

공간 부분 집합

st_disjoint는 sf 패키지에 포함된 함수로, 스페이스 객체끼리 겹치지 않는 부분을 반환하는 함수.

따라서, canterbury와 공간적으로 겹치지 않는 부분에 대해서만 nz_height 객체에서 값을 선택

canterbury  <-  nz %>% filter(Name == "Canterbury")
canterbury_height <-  nz_height[canterbury, ]

nz_height[canterbury, 2, op = st_disjoint]
Simple feature collection with 31 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1204143 ymin: 5048309 xmax: 1822492 ymax: 5650492
Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
First 10 features:
   elevation                geometry
1       2723 POINT (1204143 5049971)
2       2820 POINT (1234725 5048309)
3       2830 POINT (1235915 5048745)
4       3033 POINT (1259702 5076570)
12      2759 POINT (1373264 5175442)
25      2756 POINT (1374183 5177165)
26      2800 POINT (1374469 5176966)
27      2788 POINT (1375422 5177253)
46      2782 POINT (1383006 5181085)
47      2905 POINT (1383486 5181270)
plot(nz_height[canterbury, 2, op = st_disjoint])

st_intersects() 를 활용한 공간 부분 집합 추출

sel_sgbp <-  st_intersects(x = nz_height, y = canterbury)

sel_logical  <-  lengths(sel_sgbp) > 0

canterbury_height2 <-  nz_height[sel_logical, ]

canterbury_height3  <-  nz_height %>%
  filter(st_intersects(x = ., y = canterbury, sparse = FALSE))
Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
ℹ Please use one dimensional logical vectors instead.
ℹ The deprecated feature was likely used in the dplyr package.
  Please report the issue at <https://github.com/tidyverse/dplyr/issues>.
class(sel_sgbp)
[1] "sgbp" "list"
st_intersects(x = nz_height, y = canterbury)
Sparse geometry binary predicate list of length 101, where the
predicate was `intersects'
first 10 elements:
 1: (empty)
 2: (empty)
 3: (empty)
 4: (empty)
 5: 1
 6: 1
 7: 1
 8: 1
 9: 1
 10: 1

위상관계

# create a polygon
a_poly  <-  st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a <-   st_sfc(a_poly)
# create a line
l_line <-   st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2))
l <-   st_sfc(l_line)
# create points
p_matrix <- matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi <- st_multipoint(x = p_matrix)
p <- st_cast(st_sfc(p_multi), "POINT")

plot(a, col = c("gray"), border = c("red"))
plot(l,add = T)
plot(p,add = T)
box(col="black")

axis(side = 1, at = seq(-1.0, 1.0, 0.5), tck = 0.02)
axis(side = 2, at = seq(-1, 1, 0.5), tck = 0.02, las=1)
text(p_matrix,pos=1)

polygon과 point의 겹치는 부분을 반환

st_intersects() 겹치는 부분

st_intersects(p, a)
Sparse geometry binary predicate list of length 4, where the predicate
was `intersects'
 1: 1
 2: 1
 3: (empty)
 4: (empty)
st_intersects(p, a,sparse = F)[,1]
[1]  TRUE  TRUE FALSE FALSE

st_within() 완전히 위에 있는 부분

st_within(p, a, sparse = FALSE)[, 1]
[1]  TRUE FALSE FALSE FALSE

st_touches() 테두리만 반환

st_touches(p, a, sparse = FALSE)[, 1]
[1] FALSE  TRUE FALSE FALSE

st_is_within_distance() 는 삼각형에서 주어진 거리보다 가까운 객체들을 반환

st_is_within_distance(p, a, dist = 0.9,sparse = F)
      [,1]
[1,]  TRUE
[2,]  TRUE
[3,] FALSE
[4,]  TRUE

Spartial data operations

OSM : OpenStreetMap

cycle_hire과 cycle_hire_osm은 겹치는 포인트가 없음

plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")

any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
[1] FALSE
tmap_mode("view")
tmap mode set to interactive viewing
tm_basemap("Stamen.Terrain") +
  tm_shape(cycle_hire) + 
  tm_symbols(col = "red", shape = 16, size = 0.5, alpha = .5) + 
  tm_shape(cycle_hire_osm) +
  tm_symbols(col = "blue", shape = 16, size = 0.5, alpha = .5) + 
  tm_tiles("Stamen.TonerLabels")
tmap_mode("plot")
tmap mode set to plotting

st_is_within_distance()

위경도가 아닌 투영좌표계로 변경

cycle_hire_P <-  st_transform(cycle_hire, 27700)
cycle_hire_osm_P <-  st_transform(cycle_hire_osm, 27700)
sel <-  st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20)
summary(lengths(sel) > 0)
   Mode   FALSE    TRUE 
logical     304     438 
sum(lengths(sel) > 0)
[1] 438
mean(lengths(sel) > 0)
[1] 0.5902965
z = st_join(cycle_hire_P, cycle_hire_osm_P,
            join = st_is_within_distance, dist = 20)

nrow(cycle_hire)
[1] 742
nrow(z)
[1] 762
z = z %>% 
  group_by(id) %>% 
  summarize(capacity = mean(capacity))

nrow(z)
[1] 742
plot(cycle_hire_osm["capacity"])

plot(z["capacity"])

0508

aggregate() : group by spatial data

nz_height
Simple feature collection with 101 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1204143 ymin: 5048309 xmax: 1822492 ymax: 5650492
Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
First 10 features:
   t50_fid elevation                geometry
1  2353944      2723 POINT (1204143 5049971)
2  2354404      2820 POINT (1234725 5048309)
3  2354405      2830 POINT (1235915 5048745)
4  2369113      3033 POINT (1259702 5076570)
5  2362630      2749 POINT (1378170 5158491)
6  2362814      2822 POINT (1389460 5168749)
7  2362817      2778 POINT (1390166 5169466)
8  2363991      3004 POINT (1372357 5172729)
9  2363993      3114 POINT (1372062 5173236)
10 2363994      2882 POINT (1372810 5173419)
nz_avheight  <-  aggregate(x = nz_height, by = nz, FUN = mean)
plot(nz_avheight[2])

#group_by() 와 summarize()함수 사용
nz_avheight2  <-  nz %>%
  st_join(nz_height) %>%
  group_by(Name) %>%
  summarize(elevation = mean(elevation, na.rm = TRUE))
plot(nz_avheight2[2])

Simplification

st_simplify

seine_simp  <-  st_simplify(seine, dTolerance = 2000)  # 2000 m의 허용오차 값으로 단순화
plot(seine)

plot(seine_simp)

object.size(seine)
18096 bytes
object.size(seine_simp)
9112 bytes

differnce between st_simplify(), ms_simplify()

us_states
Simple feature collection with 49 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -124.7042 ymin: 24.55868 xmax: -66.9824 ymax: 49.38436
Geodetic CRS:  NAD83
First 10 features:
   GEOID        NAME   REGION             AREA total_pop_10 total_pop_15
1     01     Alabama    South 133709.27 [km^2]      4712651      4830620
2     04     Arizona     West 295281.25 [km^2]      6246816      6641928
3     08    Colorado     West 269573.06 [km^2]      4887061      5278906
4     09 Connecticut Norteast  12976.59 [km^2]      3545837      3593222
5     12     Florida    South 151052.01 [km^2]     18511620     19645772
6     13     Georgia    South 152725.21 [km^2]      9468815     10006693
7     16       Idaho     West 216512.66 [km^2]      1526797      1616547
8     18     Indiana  Midwest  93648.40 [km^2]      6417398      6568645
9     20      Kansas  Midwest 213037.08 [km^2]      2809329      2892987
10    22   Louisiana    South 122345.76 [km^2]      4429940      4625253
                         geometry
1  MULTIPOLYGON (((-88.20006 3...
2  MULTIPOLYGON (((-114.7196 3...
3  MULTIPOLYGON (((-109.0501 4...
4  MULTIPOLYGON (((-73.48731 4...
5  MULTIPOLYGON (((-81.81169 2...
6  MULTIPOLYGON (((-85.60516 3...
7  MULTIPOLYGON (((-116.916 45...
8  MULTIPOLYGON (((-87.52404 4...
9  MULTIPOLYGON (((-102.0517 4...
10 MULTIPOLYGON (((-92.01783 2...
us_states2163 <-  st_transform(us_states, 2163)
us_states2163
Simple feature collection with 49 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2036903 ymin: -2112380 xmax: 2521545 ymax: 731357.6
Projected CRS: NAD27 / US National Atlas Equal Area
First 10 features:
   GEOID        NAME   REGION             AREA total_pop_10 total_pop_15
1     01     Alabama    South 133709.27 [km^2]      4712651      4830620
2     04     Arizona     West 295281.25 [km^2]      6246816      6641928
3     08    Colorado     West 269573.06 [km^2]      4887061      5278906
4     09 Connecticut Norteast  12976.59 [km^2]      3545837      3593222
5     12     Florida    South 151052.01 [km^2]     18511620     19645772
6     13     Georgia    South 152725.21 [km^2]      9468815     10006693
7     16       Idaho     West 216512.66 [km^2]      1526797      1616547
8     18     Indiana  Midwest  93648.40 [km^2]      6417398      6568645
9     20      Kansas  Midwest 213037.08 [km^2]      2809329      2892987
10    22   Louisiana    South 122345.76 [km^2]      4429940      4625253
                         geometry
1  MULTIPOLYGON (((1076896 -10...
2  MULTIPOLYGON (((-1379163 -1...
3  MULTIPOLYGON (((-759892.8 -...
4  MULTIPOLYGON (((2148004 241...
5  MULTIPOLYGON (((1855614 -20...
6  MULTIPOLYGON (((1311335 -99...
7  MULTIPOLYGON (((-1298505 24...
8  MULTIPOLYGON (((1033794 -28...
9  MULTIPOLYGON (((-175289.2 -...
10 MULTIPOLYGON (((778807 -166...
us_states_simp1 <-  st_simplify(us_states2163, dTolerance = 100000)
plot(us_states)

plot(us_states_simp1)

# proportion of points to retain (0-1; default 0.05)

us_states2163$AREA <-  as.numeric(us_states2163$AREA)
us_states_simp2 <-  rmapshaper::ms_simplify(us_states2163, keep = 0.01,
                                          keep_shapes = TRUE)
plot(us_states_simp2)

Centroids

st_centroid()

지리적 개체의 중심을 식별, 때떄로 지리적 중심이 상위 개체의 경계를 벗어나는 경우가 발생(검은점)

st_point_on_surface()

상위 개체 위에 중심이 생성(빨간점)

nz_centroid  <-  st_centroid(nz)
Warning: st_centroid assumes attributes are constant over geometries
nz_pos <- st_point_on_surface(nz)
Warning: st_point_on_surface assumes attributes are constant over geometries
plot(st_geometry(nz))
plot(nz_centroid ,add=T, col="black")
Warning in plot.sf(nz_centroid, add = T, col = "black"): ignoring all but the
first attribute
plot(nz_pos ,add=T, col="red")
Warning in plot.sf(nz_pos, add = T, col = "red"): ignoring all but the first
attribute

seine_centroid <-  st_centroid(seine)
Warning: st_centroid assumes attributes are constant over geometries
seine_pos <- st_point_on_surface(seine)
Warning: st_point_on_surface assumes attributes are constant over geometries
plot(st_geometry(seine))
plot(seine_centroid ,add=T, col="black")
plot(seine_pos ,add=T, col="red")

Buffers

버퍼: 기하학적 특징의 주어진 거리 내 영역을 나타내는 다각형(입력이 점, 선 또는 다각형인지 상관없이 다각형으로 출력)

seine_buff_5km  <-  st_buffer(seine, joinStyle = "ROUND", dist = 5000)
seine_buff_50km <-  st_buffer(seine, dist = 20000)


plot(seine,col="black")
plot(seine_buff_5km,col = "red",add=T)

plot(seine,col="black")
plot(seine_buff_50km,col = "red",add=T)

Affine transformations

왜곡되거나 잘못 투영된 지도를 기반으로 생성된 지오메트리를 재투영하거나 개선할 때 많은 Affine 변환이 적용

이동하면 맵 단위로 모든 포인트가 동일한 거리만큼 이동 모든 y 좌표를 북쪽으로 100,000미터 이동하지만 x 좌표는 그대로 유지

nz_sfc  <-  st_geometry(nz)
nz_shift <-  nz_sfc + c(0, 100000)

plot(nz_sfc)
plot(nz_shift,add=T,col = "red")

배율 조정

개체를 요소만큼 확대하거나 축소

모든 기하 도형의 토폴로지 관계를 그대로 유지하면서 원점 좌표와 관련된 모든 좌표 값을 늘리거나 줄일수 있음

중심점을 기준으로 기하 도형의 차이 만큼을 늘리고 줄인 다음 다시 중심점을 더해줌

nz_centroid_sfc <-  st_centroid(nz_sfc)
nz_scale  <-  (nz_sfc - nz_centroid_sfc) * 0.5 + nz_centroid_sfc

plot(nz_sfc)
plot(nz_scale,col = "red", add=T)

회전 : 2차원 좌표의 회전하기 위한

회전변환 행렬

반시계 방향으로 회전하는 회전변환 행렬

rotation  <-  function(a){
  r = a * pi / 180 #degrees to radians
  matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
} 

nz_rotate <-  (nz_sfc - nz_centroid_sfc) * rotation(30) + nz_centroid_sfc

plot(nz_sfc)
plot(nz_rotate, add=T, col="red")

Clipping

  • 공간 클리핑은 영향을 받는 일부 형상의 지오메트리 열의 변경을 수반하는 공간 부분 집합의 한 형태

  • 클리핑은 선, 다각형 및 해당 ‘다중’ 등점보다 복잡한 형상에만 적용

  • 두 개의 겹치는 원은 중심점이 서로 1 만큼 떨어져 있고 반지름은 1인 원

b  <-  st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points
b <-  st_buffer(b, dist = 1) # convert points to circles
plot(b, border = "grey")
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3) # add text

st_intersection(), st_difference(x,y), st_difference(y,x), st_union(), st_sym_difference()

x  <-  b[1]
y <-  b[2]

x_and_y <-  st_intersection(x, y)
plot(b, border = "grey")
plot(x_and_y, col = "lightgrey", border = "grey", add = TRUE)

x_dif_y <-  st_difference(x,y)
plot(b, border = "grey")
plot(x_dif_y, col = "lightgrey", border = "grey", add = TRUE) 

y_dif_x <-  st_difference(y,x)
plot(b, border = "grey")
plot(y_dif_x, col = "lightgrey", border = "grey", add = TRUE) 

x_union_y <-  st_union(x,y)
plot(b, border = "grey")
plot(x_union_y, col = "lightgrey", border = "grey", add = TRUE)

x_sdif_y <-  st_sym_difference(x,y)
plot(b, border = "grey")
plot(x_sdif_y, col = "lightgrey", border = "grey", add = TRUE)

Subsetting(부분집합) and clipping

  • 클리핑 오브젝트는 지오메트리를 변경할 수 있지만 오브젝트의 부분 집합을 지정할 수도 있으며 클리핑/하위 설정 오브젝트와 교차하는 피쳐만 반환할 수도 있음

  • 이 점을 설명하기 위해 아래 그림은 원 x와 y의 경계 상자를 덮는 점들을 부분집합

  • 어떤 점은 원 하나 안에 있고, 어떤 점은 둘 다 안에 있고, 어떤 점은 둘 다 안에 있음

  • st_sample() : x와 y의 범위 내에서 점들의 간단한 무작위 분포를 생성하기 위해 아래에 사용

  • 그 결과 아래 그램에 표시된 출력이 나타나며, x와 y 둘 다와 교차하는 점만을 반환하기 위해 점들을 부분집합하는 방법은 무엇인가?

bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
x_and_y = st_intersection(x, y)

plot(b, border = "grey")
plot(p,add=T)

#1번째방법
p_xy1 <-  p[x_and_y]
plot(b, border = "grey")
plot(p,add=T)
plot(p_xy1,add=T,col="red")

#2번째방법
p_xy2 <-  st_intersection(p, x_and_y)
plot(p_xy2,add=T,col="blue")

#3번째방법
sel_p_xy <-  st_intersects(p, x, sparse = FALSE)[, 1] &
  st_intersects(p, y, sparse = FALSE)[, 1]
p_xy3 <-  p[sel_p_xy]
plot(p_xy3,add=T,col="green")

Geometry unions

aggregate()

regions  <-  aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
                    FUN = sum, na.rm = TRUE)
plot(regions[2])

us_states %>% View()

groupby() & summarize()

regions2  <-  us_states %>% 
  group_by(REGION) %>%
  summarize(pop = sum(total_pop_15, na.rm = TRUE))
plot(regions2[2])

st_union() 공간정보만 가져오기

us_west  <-  us_states[us_states$REGION == "West", ]
us_west_union  <-  st_union(us_west)
plot(us_west_union)

texas <-  us_states[us_states$NAME == "Texas", ]
texas_union  <-  st_union(us_west_union, texas)
plot(texas_union)

Type transformation(유형변환)

st_cast() 로 구현

multipoint <-  st_multipoint(matrix(c(1, 3, 5, 
                                      1, 3, 1), ncol = 2))
linestring <-  st_cast(multipoint, "LINESTRING")
polyg <-  st_cast(multipoint, "POLYGON")

plot(multipoint)

plot(linestring)

plot(polyg,col = "lightblue")

st_length(linestring)    #길이계산
[1] 5.656854
st_area(polyg)           #면적계산
[1] 4

st_cast() multilinestring to seperation of linstrings

multilinestring_list <-  list(matrix(c(1, 4, 5, 3), ncol = 2), 
                            matrix(c(4, 4, 4, 1), ncol = 2),
                            matrix(c(2, 4, 2, 2), ncol = 2))
multilinestring <- st_multilinestring((multilinestring_list))
multilinestring_sf <-  st_sf(geom = st_sfc(multilinestring))
multilinestring_sf
Simple feature collection with 1 feature and 0 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
CRS:           NA
                            geom
1 MULTILINESTRING ((1 5, 4 3)...
plot(multilinestring_sf)

linestring_sf2 = st_cast(multilinestring_sf, "LINESTRING")
linestring_sf2
Simple feature collection with 3 features and 0 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
CRS:           NA
                   geom
1 LINESTRING (1 5, 4 3)
2 LINESTRING (4 4, 4 1)
3 LINESTRING (2 2, 4 2)
linestring_sf2$name  <-  c("Riddle Rd", "Marshall Ave", "Foulke St")
linestring_sf2$length <-  st_length(linestring_sf2)
linestring_sf2
Simple feature collection with 3 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
CRS:           NA
                   geom         name   length
1 LINESTRING (1 5, 4 3)    Riddle Rd 3.605551
2 LINESTRING (4 4, 4 1) Marshall Ave 3.000000
3 LINESTRING (2 2, 4 2)    Foulke St 2.000000
plot(linestring_sf2[2])

0515 Cpt.9

Tmap package basics

tmap_mode("plot")
tmap mode set to plotting
tm_shape(nz)+
  tm_fill()

tm_shape(nz)+
  tm_borders()

tm_shape(nz)+
  tm_borders()+
  tm_fill()

Map objects

map_nz <- tm_shape(nz) + tm_polygons()

map_nz1 <-  map_nz +
  tm_shape(nz_elev) + 
  tm_raster(alpha = 0.7)

nz_water <- st_union(nz) %>% 
  st_buffer(22200) %>% 
  st_cast(to = "LINESTRING")

map_nz2 <- map_nz1 +
  tm_shape(nz_water) + 
  tm_lines()

map_nz3 <-  map_nz2 +
  tm_shape(nz_height) + 
  tm_dots()

tmap_arrange(map_nz1, map_nz2, map_nz3)
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
Some legend labels were too wide. These labels have been resized to 0.62, 0.55, 0.55, 0.55, 0.55, 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
Some legend labels were too wide. These labels have been resized to 0.62, 0.55, 0.55, 0.55, 0.55, 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
Some legend labels were too wide. These labels have been resized to 0.62, 0.55, 0.55, 0.55, 0.55, 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Aesthetics

ma1 <-  tm_shape(nz) + tm_fill(col = "red")
ma2 <-  tm_shape(nz) + tm_fill(col = "red", alpha = 0.3)
ma3 <-  tm_shape(nz) + tm_borders(col = "blue")
ma4 <-  tm_shape(nz) + tm_borders(lwd = 3)
ma5 <-  tm_shape(nz) + tm_borders(lty = 2)
ma6 <-  tm_shape(nz) + tm_fill(col = "red", alpha = 0.3) +
  tm_borders(col = "blue", lwd = 3, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)

plot(st_geometry(nz), col = nz$Land_area)  # works

tm_shape(nz) + tm_fill(col = "Land_area")

Set legend title

map_nza  <-  tm_shape(nz) +
  tm_fill(col = "Land_area", title = expression("Area (km"^2*")")) + 
  tm_borders()

map_nza

Set Legend breaks & Palette

a <- tm_shape(nz) + tm_polygons(col = "Median_income")
b <- tm_shape(nz) + tm_polygons(col = "Median_income", breaks = c(0, 3, 4, 5) * 10000)
c <- tm_shape(nz) + tm_polygons(col = "Median_income", n = 10)
d <- tm_shape(nz) + tm_polygons(col = "Median_income", palette = "BuGn")

tmap_arrange(a,b,c,d)

  • style = "pretty" :기본 설정은 가능한 경우 정수로 반올림하고 간격을 균등하게 유지.

  • style = "equal" : 입력 값을 동일한 범위의 빈으로 나누고 균일한 분포의 변수에 적합(결과 맵이 색상 다양성이 거의 없을 수 있으므로 분포가 치우친 변수에는 권장하지 않음).

  • style = "quantile" : 동일한 수의 관찰이 각 범주에 포함되도록 함(빈 범위가 크게 다를 수 있다는 잠재적인 단점이 있음).

  • style = "jenks" : 데이터에서 유사한 값의 그룹을 식별하고 범주 간의 차이를 최대화

  • style = "cont" : 연속 색상 필드에 많은 색상을 표시하고 연속 래스터에 특히 적합

  • style = "cat" : 범주 값을 나타내도록 설계되었으며 각 범주가 고유한 색상을 받도록 함

a <- tm_shape(nz) + tm_polygons(col = "Median_income", style = "pretty")
b <- tm_shape(nz) + tm_polygons(col = "Median_income", style = "equal")
c <- tm_shape(nz) + tm_polygons(col = "Median_income", style = "quantile")
d <- tm_shape(nz) + tm_polygons(col = "Median_income", style = "jenks")

tmap_arrange(a,b,c,d)

tm_shape(nz) + tm_polygons("Population", palette = "Blues")
Some legend labels were too wide. These labels have been resized to 0.64, 0.59, 0.59. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

tm_shape(nz) + tm_polygons("Population", palette = "YlOrBr")
Some legend labels were too wide. These labels have been resized to 0.64, 0.59, 0.59. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Layouts

a <- map_nz + 
  tm_compass(type = "8star", position = c("left", "top")) +
  tm_scale_bar(breaks = c(0, 100, 200), text.size = 1)

b <- map_nz + tm_layout(title = "New Zealand")
c <- map_nz + tm_layout(scale = 5)
d <- map_nz + tm_layout(bg.color = "lightblue")
e <- map_nz + tm_layout(frame = FALSE)

tmap_arrange(a,b,c,d,e,nrow = 1)

tm_layout() 의 다양한 옵션

  • Frame width (frame.lwd) : 프레임 너비

  • an option to allow double lines (frame.double.line ) : 이중선 허용 옵션

  • Margin settings including outer.margin and inner.margin : 여백 설정

  • Font settings controlled by fontface and fontfamily : 글꼴 설정

  • legend.show (whether or not to show the legend) : 범례 표시 여부

  • multiple choice settings such as legend.position :범례 위치 변경

tm_style()

a <- map_nza + tm_style("bw")
b <- map_nza + tm_style("classic")
c <- map_nza + tm_style("cobalt")
d <- map_nza + tm_style("col_blind")

tmap_arrange(a,b,c,d)

Faceted maps

  • Faceted maps은 나란히 배열된 많은 맵으로 구성

  • Faceted maps을 사용하면 시간과 같은 다른 변수와 관련하여 공간 관계가 어떻게 변하는지 시각화할 수 있음

urb_1970_2030  <-  urban_agglomerations %>% 
  filter(year %in% c(1970, 1990, 2010, 2030))

tm_shape(world) +
  tm_polygons() +
  tm_shape(urb_1970_2030) +
  tm_symbols(col = "black", border.col = "white", size = "population_millions") +
  tm_facets(by = "year", nrow = 2, free.coords = FALSE)

urb_anim  <-  tm_shape(world) + tm_polygons() + 
  tm_shape(urban_agglomerations) + tm_dots(size = "population_millions") +
  tm_facets(along = "year", free.coords = FALSE)

tmap_animation(urb_anim, filename = "urb_anim.gif", delay = 25)
Creating frames
=========
====
=====
====
=====
====
=====
====
====
=====
====
=====
====
=====
====
=====
====

Creating animation
Animation saved to C:\sanai_sungil\basics_example\urb_anim.gif

Inset maps

  • Inset maps는 기본 지도 내부 또는 옆에 렌더링되는 더 작은 지도

  • 첫 번째 단계는 관심 영역을 정의하는 것인데, 이는 새로운 공간 객체를 생성하여 수행할 수 있음 nz_region.

nz_region <-  st_bbox(c(xmin = 1340000, xmax = 1450000,
                      ymin = 5130000, ymax = 5210000),
                    crs = st_crs(nz_height)) %>%  st_as_sfc()
# 뉴질랜드의 서던 알프스 지역을 보여주는 기본 지도를 만듬
nz_height_map  <-  tm_shape(nz_elev, bbox = nz_region) +
  tm_raster(style = "cont", palette = "YlGn", legend.show = TRUE) +
  tm_shape(nz_height) + 
  tm_symbols(shape = 2, col = "red", size = 1) +
  tm_scale_bar(position = c("left", "bottom"))

# 삽입 맵 생성으로 구성
nz_map <- tm_shape(nz) + 
  tm_polygons() +
  tm_shape(nz_height) + 
  tm_symbols(shape = 2, col = "red", size = 0.1) + 
  tm_shape(nz_region) + 
  tm_borders(lwd = 3)

#  `viewport()` 패키지의 함수를 사용하여 두 맵을 결합
nz_height_map
stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
print(nz_map, vp = viewport(0.8, 0.27, width = 0.5, height = 0.5))

Interactive maps

  • 대화형 지도는 데이터 세트를 새로운 차원으로 끌어올릴 수 있음

  • ’웹 지도’에 오버레이된 지리 데이터 세트의 모든 부분을 이동하고 확대하는 기능

  • 지도를 기울이고 회전하는 기능과 사용자가 이동 및 확대/축소할 때 자동으로 업데이트됨

  • tmap, mapview,mapdeck, leaflet 을 가지고 대화형 지도 표현 가능

tmap_mode("view") #interactive mode가 on상태임
tmap mode set to interactive viewing
map_nz <- tm_shape(nz)+
  tm_polygons()
  
map_nz

mapdeck

library(mapdeck)

Attaching package: 'mapdeck'
The following object is masked from 'package:tibble':

    add_column
set_token("pk.eyJ1Ijoic3VuZ2lsZW8iLCJhIjoiY2xoYTRwbXEzMGR6eTNkbXBoZnluNXdyYSJ9.Id1fKIbhtvA9Mrnyo_1JQA")
crash_data = read.csv("https://git.io/geocompr-mapdeck")
crash_data = na.omit(crash_data)

ms = mapdeck_style("dark")

mapdeck(style = ms, 
        pitch = 45, 
        location = c(0, 52), 
        zoom = 4) %>%
  add_grid(data = crash_data, 
           lat = "lat", 
           lon = "lng", 
           cell_size = 1000,
           elevation_scale = 50, 
           layer_id = "grid_layer",
           colour_range = viridisLite::plasma(6))
Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite

more interactive map examples

Geospatial Visualization